/*

  Simple cross-product, outputs in correct triangular form.

  - Coalescing memory access in all reads
  - No memory coalscing in writes
  - Shared memory reads of type float2 to reduce global memory traffic
  - Each thread works on a 2x1 square of data

 */

// Read in column in first warp as float2, row in second warp
#define LOAD_TWO_BY_ONE(s)						\
  {*(input##s##_p) = *array_p;						\
    CUBE_ADD_BYTES(4);}

// read in shared data as individual floats to avoid bank conflicts

#define TWO_BY_ONE_COMPUTE(s)						\
  {float col1Xreal = input[s][colIdx];					\
  float col1Ximag = input[s][colIdx + 1];			        \
  float col1Yreal = input[s][colIdx + 2];				\
  float col1Yimag = input[s][colIdx + 3];				\
  float row1Xreal = input[s][rowIdx];					\
  float row1Ximag = input[s][rowIdx + 1];				\
  float row1Yreal = input[s][rowIdx + 2];				\
  float row1Yimag = input[s][rowIdx + 3];				\
  float row2Xreal = input[s][rowIdx + 4];				\
  float row2Ximag = input[s][rowIdx + 5];				\
  float row2Yreal = input[s][rowIdx + 6];				\
  float row2Yimag = input[s][rowIdx + 7];				\
  sum11XXreal += row1Xreal * col1Xreal;					\
  sum11XXreal += row1Ximag * col1Ximag;					\
  sum11XXimag += row1Ximag * col1Xreal;					\
  sum11XXimag -= row1Xreal * col1Ximag;					\
  sum11XYreal += row1Xreal * col1Yreal;					\
  sum11XYreal += row1Ximag * col1Yimag;					\
  sum11XYimag += row1Ximag * col1Yreal;					\
  sum11XYimag -= row1Xreal * col1Yimag;					\
  sum11YXreal += row1Yreal * col1Xreal;					\
  sum11YXreal += row1Yimag * col1Ximag;					\
  sum11YXimag += row1Yimag * col1Xreal;					\
  sum11YXimag -= row1Yreal * col1Ximag;					\
  sum11YYreal += row1Yreal * col1Yreal;					\
  sum11YYreal += row1Yimag * col1Yimag;					\
  sum11YYimag += row1Yimag * col1Yreal;					\
  sum11YYimag -= row1Yreal * col1Yimag;					\
  sum21XXreal += row2Xreal * col1Xreal;					\
  sum21XXreal += row2Ximag * col1Ximag;					\
  sum21XXimag += row2Ximag * col1Xreal;					\
  sum21XXimag -= row2Xreal * col1Ximag;					\
  sum21XYreal += row2Xreal * col1Yreal;					\
  sum21XYreal += row2Ximag * col1Yimag;					\
  sum21XYimag += row2Ximag * col1Yreal;					\
  sum21XYimag -= row2Xreal * col1Yimag;					\
  sum21YXreal += row2Yreal * col1Xreal;					\
  sum21YXreal += row2Yimag * col1Ximag;					\
  sum21YXimag += row2Yimag * col1Xreal;					\
  sum21YXimag -= row2Yreal * col1Ximag;					\
  sum21YYreal += row2Yreal * col1Yreal;					\
  sum21YYreal += row2Yimag * col1Yimag;					\
  sum21YYimag += row2Yimag * col1Yreal;					\
  sum21YYimag -= row2Yreal * col1Yimag;}

/*  sum11XXreal += row1Xreal * col1Xreal + row1Ximag * col1Ximag;	\
  sum11XXimag += row1Ximag * col1Xreal - row1Xreal * col1Ximag;		\
  sum11XYreal += row1Xreal * col1Yreal + row1Ximag * col1Yimag;		\
  sum11XYimag += row1Ximag * col1Yreal - row1Xreal * col1Yimag;		\
  sum11YXreal += row1Yreal * col1Xreal + row1Yimag * col1Ximag;		\
  sum11YXimag += row1Yimag * col1Xreal - row1Yreal * col1Ximag;		\
  sum11YYreal += row1Yreal * col1Yreal + row1Yimag * col1Yimag;		\
  sum11YYimag += row1Yimag * col1Yreal - row1Yreal * col1Yimag;		\
  sum21XXreal += row2Xreal * col1Xreal + row2Ximag * col1Ximag;		\
  sum21XXimag += row2Ximag * col1Xreal - row2Xreal * col1Ximag;		\
  sum21XYreal += row2Xreal * col1Yreal + row2Ximag * col1Yimag;		\
  sum21XYimag += row2Ximag * col1Yreal - row2Xreal * col1Yimag;		\
  sum21YXreal += row2Yreal * col1Xreal + row2Yimag * col1Ximag;		\
  sum21YXimag += row2Yimag * col1Xreal - row2Yreal * col1Ximag;		\
  sum21YYreal += row2Yreal * col1Yreal + row2Yimag * col1Yimag;		\
  sum21YYimag += row2Yimag * col1Yreal - row2Yreal * col1Yimag;}*/


CUBE_KERNEL(shared2x1float, cuFloatComplex *array_d, cuFloatComplex *product_d, 
	    const int write)
{
  CUBE_START;

  //get local thread ID
  int ty = threadIdx.y;
  int tx = threadIdx.x;
  int tid = ty*TILE_WIDTH + tx;

  //set frequency number from blockIdx.y
  int f = blockIdx.y;

  int Row, Col, blockX, blockY;
  CUBE_DEVICE_CALL(findPosition, Col, Row, blockX, blockY);
  Row*=2;

  //declare shared memory for input coalescing
  #define PAD 1
  __shared__ float input[2][4*(2*TILE_WIDTH+PAD)]; // 4* for float4, 2* for 2x2 tile size

  //instantiate sum variables
  float sum11XXreal = 0.0, sum11XXimag = 0.0;
  float sum11XYreal = 0.0, sum11XYimag = 0.0;
  float sum11YXreal = 0.0, sum11YXimag = 0.0;
  float sum11YYreal = 0.0, sum11YYimag = 0.0;
  float sum21XXreal = 0.0, sum21XXimag = 0.0;
  float sum21XYreal = 0.0, sum21XYimag = 0.0;
  float sum21YXreal = 0.0, sum21YXimag = 0.0;
  float sum21YYreal = 0.0, sum21YYimag = 0.0;

  float *array_p = (float*)array_d + tid;
  if (tid < 4*TILE_WIDTH) {
    array_p += (f*NSTATION*NPOL + blockX*TILE_WIDTH*NPOL)*2; //we're dealing with floats instead of float2s
  } else {
    array_p += (f*NSTATION*NPOL + blockY*TILE_WIDTH*NPOL)*2 - 4*TILE_WIDTH;    
  }

  // introduce colIdx and rowIdx variables for conflict-free reading
  int colIdx = (tx/8)*(2*TILE_WIDTH+PAD) + 4*(tx%8);
  int rowIdx = (ty/4)*(2*TILE_WIDTH+PAD) + 8*(ty%4) + 2*(2*TILE_WIDTH + PAD);

  // introduce new tid value for conflict-free writing
  tid = (2*TILE_WIDTH + PAD)*(ty/2) + TILE_WIDTH*(ty%2) + tx;

  float *input0_p = input[0] + tid;
  float *input1_p = input[1] + tid;

  LOAD_TWO_BY_ONE(0);

  //#pragma unroll 2
  for (int t=0; t<NTIME-2; t+=2){
    __syncthreads();
    TWO_BY_ONE_COMPUTE(0);
    array_p += NFREQUENCY*NSTATION*NPOL*2;

    LOAD_TWO_BY_ONE(1);
    __syncthreads();
    TWO_BY_ONE_COMPUTE(1);
    array_p += NFREQUENCY*NSTATION*NPOL*2;

    LOAD_TWO_BY_ONE(0);
  }

  __syncthreads();

  TWO_BY_ONE_COMPUTE(0);
  array_p += NFREQUENCY*NSTATION*NPOL*2;

  LOAD_TWO_BY_ONE(1);
  __syncthreads();
  TWO_BY_ONE_COMPUTE(1);

  if (write)
    CUBE_DEVICE_CALL(write2x1, Col, Row, product_d, sum11XXreal, sum11XXimag, sum11XYreal, 
		     sum11XYimag, sum11YXreal, sum11YXimag, sum11YYreal, sum11YYimag, 
		     sum21XXreal, sum21XXimag, sum21XYreal, sum21XYimag, sum21YXreal, 
		     sum21YXimag, sum21YYreal, sum21YYimag);

  if (Col < Row+1) {
    CUBE_ADD_BYTES(64);
    CUBE_ADD_FLOPS(NTIME*64);
  } else {
    CUBE_ADD_BYTES(32);
    CUBE_ADD_FLOPS(NTIME*32);
  }

  CUBE_END;
}

